// Molecule.cpp: implementation of the Molecule class.
//
//////////////////////////////////////////////////////////////////////

#include "Molecule.h"
#include "nr_c/nr.h"
#include <vector>

#include <iomanip>
using namespace std;
//#include <cstdio>

#define PI   3.1415926535897932385

#define CONV 0.01745329252

// If running famotodine or capsaicin, we must
// set this TRANSLATION to be true. If running
// chlorothiazide must set this TRANSLATION to be
// false.
#define TRANSLATION true

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////



Molecule::Molecule(int i_NInput, char *ZMatrixFile)
{

  NInput = i_NInput;

  newA(InputValue, NInput);
  newA(ZmatrixInputValue, NInput);
  newA(Derivative, NInput);


  SAFEifstream infile(ZMatrixFile);

  	
  double len1, len2, len3, phi1, phi2, phi3;
  	
  infile >> len1 >> len2 >> len3 >> phi1 >> phi2 >> phi3;

  // abc and abg are used in PrintCartesians() method

  abc[0] = len1; abc[1] = len2; abc[2] = len3; 
  abg[0] = phi1; abg[1] = phi2; abg[2] = phi3;


  phi1 *= CONV;	phi2 *= CONV;	phi3 *= CONV;

  double   s2 = (cos(phi3) - cos(phi1)*cos(phi2))/sin(phi1);
  double   s3 = cos(phi2);
  double   s1 = sqrt(1 - s2*s2 - s3*s3);
  double   s4 = sin(phi1);
  double   s5 = cos(phi1);

  UnitCell1 = Position(s1,   s2,  s3);
  UnitCell2 = Position(0.0,  s4,  s5);
  UnitCell3 = Position(0.0, 0.0, 1.0);

  UnitCell1 *= len1;
  UnitCell2 *= len2;
  UnitCell3 *= len3;

  double UnitVolume = UnitCell1*(UnitCell2^UnitCell3);

  WaveVector1 = (UnitCell2^UnitCell3)/UnitVolume;
  WaveVector2 = (UnitCell3^UnitCell1)/UnitVolume;
  WaveVector3 = (UnitCell1^UnitCell2)/UnitVolume;
  	

  // read in number of molecules in asym and no. of 
  // each molecule have

  infile >> NMolecule;
  newA(NAtomEach, NMolecule);

  NAtom = 0;
  
  for (int i = 1; i <= NMolecule; i++) 
  {
	  infile >> NAtomEach[i];
	  NAtom += NAtomEach[i];
  }  

  newA(AtomPosition,	NAtom);
  newA(AtomBuilder,	NAtom);				
  newM(atomLabel, NAtom, 3);
  newM(InputOnX,		NInput, NAtom+1);

  // each molecule in the asym much have reference base

  newA(P1, NMolecule);
  newA(P2, NMolecule);
  newA(P3, NMolecule);
{
  for (int i = 1; i <= NMolecule; i++)
  {
	  P1[i] = Position( 1.0, 0.0, 0.0);
	  P2[i] = Position( 1.0, 1.0, 0.0);	
	  P3[i] = Position( 0.0, 1.0, 0.0);
  }
}

	{
	int      i=1;

	double   *pd1, *pd2, *pd3, d1, d2, d3, dummy;
	Pack     *Dd1, *Dd2, *Dd3;
	Position *A1,  *A2,  *A3,  *A4;
	int 	 i1, i2, i3;
	int      n1, n2, n3;

	
	if(TRANSLATION)
	{
		newA(Mover, NMolecule);	// a Mover for each molecule

		for (int j = 1; j <= NMolecule; j++)
		{
			i += 3;

			// setting up pointers to be passed onto Mover

			pd1 = &InputValue[1+(j-1)*3]; 
			pd2 = &InputValue[2+(j-1)*3];
			pd3 = &InputValue[3+(j-1)*3];

			Dd1 = &Derivative[1+(j-1)*3];
			Dd2 = &Derivative[2+(j-1)*3];
			Dd3 = &Derivative[3+(j-1)*3];

			Mover[j].Initialise(P1[j], P2[j], P3[j],  &P1[j],  &P2[j],  &P3[j], 
							  pd1, pd2, pd3, Dd1, Dd2, Dd3);
		}
	}

	int n = 1; // index for degrees of freedom (i.e. the variables in the FOM)
	int iAtomDelta = 0; // no. atoms read in so far in clumps of molecules

	for (int iMol = 1; iMol <= NMolecule; iMol++)
	{
		for (int iAtom = 1; iAtom <= NAtomEach[iMol]; iAtom++)
		{
			infile >> atomLabel[n] >> d1 >> i1 >> d2 >> i2 >> d3 >> i3 
						>> n1 >> n2 >> n3 >> dummy;

			if(i1)	{ 
        pd1 = &InputValue[i]; Dd1 = &Derivative[i++]; 
        ZmatrixInputValue[i-1] = d1*PI/180.0;
      }
			else	{ pd1 = NULL; Dd1 = NULL;}

			if(i2)	{ 
        pd2 = &InputValue[i]; Dd2 = &Derivative[i++]; 
        ZmatrixInputValue[i-1] = d2*PI/180.0;
      }
			else	{ pd2 = NULL; Dd2 = NULL;}

			if(i3)  { 
        pd3 = &InputValue[i]; Dd3 = &Derivative[i++]; 
        ZmatrixInputValue[i-1] = d3*PI/180.0;
      }
			else	{ pd3 = NULL; Dd3 = NULL;}

			switch(iAtom)
			{
				case 1:
					A1 = &P1[iMol];
					A2 = &P2[iMol];
					A3 = &P3[iMol];
					break;
				case 2:
					A1 = &P2[iMol];
					A2 = &P3[iMol];
					A3 = &AtomPosition[n1 + iAtomDelta];
					break;
				case 3:
					A1 = &P3[iMol];
					A2 = &AtomPosition[n2 + iAtomDelta];
					A3 = &AtomPosition[n1 + iAtomDelta];
					break;
				default:
					A1 = &AtomPosition[n3 + iAtomDelta];
					A2 = &AtomPosition[n2 + iAtomDelta];
					A3 = &AtomPosition[n1 + iAtomDelta];
					break;
				}

			A4 = &AtomPosition[n];

			d2*=CONV; d3*=CONV; 

			AtomBuilder[n] = new Atom(d1,d2,d3,pd1,pd2,pd3,Dd1,Dd2,Dd3,A1,A2,A3,A4);

			n++;
		}
		iAtomDelta += NAtomEach[iMol]; // add clump of atoms from molecule just read in
//		cout << i << endl;
  }

  // don't currently understand why I need two getline to read the spacegroup... but never mind
  infile.getline(SpaceGroup,15);
  infile.getline(SpaceGroup,15);

	}

infile.close();

BuildDependence();
}


void Molecule::BuildDependence(void)
{
int n,m,p;
int ZMatrixError=0;

double   *TestInput;
Position *StoredPosition, *Dummy, Delta;

newA(TestInput,		 NInput);
newA(StoredPosition, NAtom);

for(m=1;m<=NInput;m++)	   
	TestInput[m] = gasdev();
/*TestInput[1] = 0.0;
TestInput[2] = 0.0;
TestInput[3] = 0.0;
TestInput[4] = 0.0;
TestInput[5] = 0.0;
TestInput[6] = 0.0;
TestInput[7] = 0.7854;
TestInput[8] = 0.0;
TestInput[9] = 0.0;
TestInput[10] = 0.7854;
TestInput[11] = 0.0;
TestInput[12] = 0.0;*/
Dummy = Build(TestInput);


// to print out some random structures
//for(m=1;m<=NInput;m++)
//  cout << (180/PI)*InputValue[m] << endl;
//PrintCartesians(TestInput);
//exit(1);


for(n=1;n<=NAtom;n++)	   
	StoredPosition[n] = AtomPosition[n];


double Change;


for(m=1;m<=NInput;m++)
	{
	Change = gasdev();				


	TestInput[m] += Change;
	Dummy         = Build(TestInput);	
	TestInput[m] -= Change;


	
	for(p=1,n=1;n<=NAtom;n++)
		{
		if(AtomPosition[n]!=StoredPosition[n])
			{
			InputOnX[m][p++] = n;
			Delta = AtomPosition[n] - StoredPosition[n];

			if(fabs(Delta*Derivative[m].Vector1) >=1e-8)
				{
				ZMatrixError = 1;
				cerr << "Problem with z-matrix : ";
				cerr << "Atom " << n << " on "; 
				cerr << "Parameter " << m << "\n";
				}
			}
		}

	InputOnX[m][p] = -1;		// array terminator
	}


deleteA(TestInput);
deleteA(StoredPosition);


if(ZMatrixError) 
	{
	cerr << "\n\n\nCan not continue.\n\n\n";
	exit(1);
	}
}


Molecule::~Molecule()
{
int n;
	for(n=1;n<=NAtom;n++)	
		delete AtomBuilder[n];
	deleteA(AtomBuilder);
	deleteA(AtomPosition);	
	deleteA(Derivative);	
	deleteA(InputValue);
	deleteA(ZmatrixInputValue);
	deleteM(InputOnX);

	deleteA(P1); deleteA(P2); deleteA(P3);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////


Position* Molecule::Build(double *Value)
{
CopyA(InputValue, Value, NInput);

if(TRANSLATION)
	for (int i = 1; i <= NMolecule; i++)
		Mover[i].Build();

for(int n=1;n<=NAtom;n++)	AtomBuilder[n]->Build();	

return (AtomPosition);
}

Position* Molecule::Build(vector<double> &Value)
{
  //CopyA(InputValue, Value, NInput);
  CopyVA(InputValue, Value, NInput);

if(TRANSLATION)
	for (int i = 1; i <= NMolecule; i++)
		Mover[i].Build();

for(int n=1;n<=NAtom;n++)	AtomBuilder[n]->Build();	

return (AtomPosition);
}

void Molecule::CopyVA(double *InputValue, vector<double> &Input, int NInput)
{
  for(int i = 1; i <= NInput; i++)
  {
    InputValue[i] = Input[i-1];
  }
}


void Molecule::TransformDerivatives(Position *XDerivative, double *InputDerivative)
{
int n,m,p;

for(m=1;m<=NInput;m++)
	{
	InputDerivative[m] = 0.0;
	for(p=1;(n=InputOnX[m][p])!=-1;p++)
		{		
		InputDerivative[m] += XDerivative[n]
							  *( (Derivative[m].Vector1^AtomPosition[n]) 
								+ Derivative[m].Vector2);
		}	
	}
}

void Molecule::PrintCIF(double *Value, const double ChiSq, char name[])
{
  ofstream out;
  out.open( (string(name)+".cif").c_str() );

  PrintCIF(Value, ChiSq, out, string("data_") + string(name));

  out.close();
}

void Molecule::PrintCIF(double *Value, const double ChiSq, ofstream& out, string cifLabel)
{
  Build(Value);

  // some start info; taking from a .cif file with P1


  out << "data_" << cifLabel << endl;
  out << "_symmetry_space_group_name_H-M '" << SpaceGroup << "'" << endl;
  out << "_cell_length_a    " << abc[0] << endl;
  out << "_cell_length_b    " << abc[1] << endl;
  out << "_cell_length_c    " << abc[2] << endl;
  out << "_cell_angle_alpha   " << abg[0] << endl;
  out << "_cell_angle_beta    " << abg[1] << endl;
  out << "_cell_angle_gamma    " << abg[2] << endl;
  out << "_cell_formula_units_Z     1" << endl;
  out << "loop_" << endl;
  out << "    _atom_site_label" << endl;
  out << "    _atom_site_fract_x" << endl;
  out << "    _atom_site_fract_y" << endl;
  out << "    _atom_site_fract_z" << endl;
  out << "    _atom_site_type_symbol" << endl;


  double dummy1[3], dummy2[3]; // for translation
  Position trans_part_cart;

  for (int i = 0; i <= 2; i++)
  { 
	  dummy1[i] = 0.0;
	  dummy2[i] = 0.0;
  }

  int n = 1; // index for atom
  int iAtomDelta = 0; // no. atoms read in so far in clumps of molecules

  for (int iMol = 1; iMol <= NMolecule; iMol++)
  {
	  if (TRANSLATION)
	  {
		  trans_part_cart = Position(InputValue[1+(iMol-1)*3], InputValue[2+(iMol-1)*3], 
			  InputValue[3+(iMol-1)*3]);

		  dummy1[0] = trans_part_cart*WaveVector1;
		  dummy1[1] = trans_part_cart*WaveVector2;
		  dummy1[2] = trans_part_cart*WaveVector3;

		  for (int i = 0; i <= 2; i++)
		  {
			  dummy2[i] = dummy1[i];
			  while (dummy1[i] < 0.0) dummy1[i] += 1.0;
			  while (dummy1[i] >= 1.0) dummy1[i] -= 1.0;
			  dummy2[i] -= dummy1[i];
		  }
	  }
	  for (int iAtom = 1; iAtom <= NAtomEach[iMol]; iAtom++)
	  {
		  out << setw(4) << atomLabel[n] << n << " ";
		  out.precision(5);
		  out << setw(10) << (AtomPosition[n]*WaveVector1)-dummy2[0] << " ";
		  out << setw(10) << (AtomPosition[n]*WaveVector2)-dummy2[1] << " ";
		  out << setw(10) << (AtomPosition[n]*WaveVector3)-dummy2[2] << setw(4) << atomLabel[n] << endl;

		  n++;
	  }
	  iAtomDelta += NAtomEach[iMol];
  }

  out << "#END" << endl;
}



void Molecule::PrintCartesians(double *Value)
{
Build(Value);

double n1 = floor(AtomPosition[1]*WaveVector1);
double n2 = floor(AtomPosition[1]*WaveVector2);
double n3 = floor(AtomPosition[1]*WaveVector3);

Position OffSet(n1*UnitCell1 + n2*UnitCell2 + n3*UnitCell3);

ofstream out;
out.open("Cart.out");	
out << "Printing out Cartesians coordinates" << endl;

int n;

//Position dum;

for(n=1;n<=NAtom;n++)
	//dum = AtomPosition[n] - OffSet;
	out << n << "\t" << AtomPosition[n] - OffSet << "\n";		

out.close();
}

void Molecule::PrintFractional(double *Value, double ChiSq, char name[], int NAccept, int NTotal)
{
  Build(Value);
  double dummy1[3], dummy2[3]; // for translation
  Position trans_part_cart;

  for (int i = 0; i <= 2; i++)
  { 
	  dummy1[i] = 0.0;
	  dummy2[i] = 0.0;
  }


  ofstream out;
  out.open( (string(name)+".cssr").c_str() );

  // see Josuttis book for explanation of this format flag
  // I don't know have to turn this flag off. However this might be 
  // because if use an old header <iomanip.h> instead of <iomanip> somehow

  out.setf(ios::fixed); 
  out << " REFERENCE STRUCTURE = 00000   A,B,C =";

  out.precision(3);
  out << setw(8) << abc[0] << setw(8) << abc[1] << setw(8) << abc[2]  << endl;

  out << "   ALPHA,BETA,GAMMA =";
  out << setw(8) << abg[0] << setw(8) << abg[1] << setw(8) << abg[2]  << endl;

  out << setw(4) << NAtom << "   0  HMC solution" << endl;
  out.precision(1);
  out << " Final Chi-Sq = " << ChiSq << " ; No samples = " << NTotal <<
	  " and no. accepted = " << NAccept << endl; 


  int n = 1; // index for atom
  int iAtomDelta = 0; // no. atoms read in so far in clumps of molecules

  for (int iMol = 1; iMol <= NMolecule; iMol++)
  {
	  if (TRANSLATION)
	  {
		  trans_part_cart = Position(InputValue[1+(iMol-1)*3], InputValue[2+(iMol-1)*3], 
			  InputValue[3+(iMol-1)*3]);

		  dummy1[0] = trans_part_cart*WaveVector1;
		  dummy1[1] = trans_part_cart*WaveVector2;
		  dummy1[2] = trans_part_cart*WaveVector3;

		  for (int i = 0; i <= 2; i++)
		  {
			  dummy2[i] = dummy1[i];
			  while (dummy1[i] < 0.0) dummy1[i] += 1.0;
			  while (dummy1[i] >= 1.0) dummy1[i] -= 1.0;
			  dummy2[i] -= dummy1[i];
		  }
	  }
	  for (int iAtom = 1; iAtom <= NAtomEach[iMol]; iAtom++)
	  {

		  out << setw(4) << n << " ";
		  out.setf(ios::left);
		  out << setw(4) << atomLabel[n];
		  out.unsetf(ios::left);
		  out.precision(5);
		  out << setw(10) << (AtomPosition[n]*WaveVector1)-dummy2[0];
		  out << setw(10) << (AtomPosition[n]*WaveVector2)-dummy2[1];
		  out << setw(10) << (AtomPosition[n]*WaveVector3)-dummy2[2] << endl;

		  n++;
	  }
	  iAtomDelta += NAtomEach[iMol];
  }

  out.close();


  // Also just print out internal torsion angle to test internal torsion
  // angle distance messure

/*
  out.open((string(name) + "_internal_dof.txt").c_str());	

  for(int i = 1; i <= NInput; i++)
  {
	  if (TRANSLATION && i <= 3)
	  {
      out << Value[i] << endl;
    }
    else
    {
      out << (180/PI)*fmod(Value[i], 2.0*PI) << endl;
    }
  }
  out.close();
*/
}


void Molecule::createRandomSol(double *TestInput)
{
	int m;

	if (TRANSLATION)
	{
		TestInput[1] = ran1()*abc[0];
		TestInput[2] = ran1()*abc[1];
		TestInput[3] = ran1()*abc[2];

		for(m=4;m<=NInput;m++)	   
		  TestInput[m] = ran1()*2.0*PI;
	}
	else
	{
		for(m=1;m<=NInput;m++)	   
		  TestInput[m] = ran1()*2.0*PI;
	}
}

// move just i'th coordinate, where i=0 here refer to first 
// coordinate
void Molecule::customModifySol(double *TestInput, int i)
{
  if (i < 0)
  {
    cout << "require i >= 0 in molecule::customModifySol\n";
    return;
  }

  i = i % NInput;

  double devideBy = 3.0;


	if (TRANSLATION)
	{
    if (i < 3)
    {
		  TestInput[i+1] += abc[i] / devideBy;
    }
    else 
    {
		  TestInput[i+1] += 2.0*PI / devideBy;
    }
	}
	else
	{
		TestInput[i+1] += 2.0*PI / devideBy;
	}
}

// distance from input z-matrix values to local solution
double Molecule::distanceFromInputZmatrix(double *solution)
{
  double toReturn = 0.0;

	if (TRANSLATION)
	{
    if ( NInput <= 6 )
      return 10000.0;

    // only compare between internal torsion angles, since rotation of molecule ga ga
    // in z matrix file

    for(int i = 7; i <= NInput; i++)
    {
      double val = getDistanceBetweenAngles(solution[i], ZmatrixInputValue[i]);
      toReturn += val*val;
    }
  }
  else
  {
    if ( NInput <= 3 )
      return 20000.0;

    // only compare between internal torsion angles, since rotation of molecule ga ga
    // in z matrix file

    for(int i = 4; i <= NInput; i++)
    {
      double val = getDistanceBetweenAngles(solution[i], ZmatrixInputValue[i]);
      toReturn += val*val;
    }
  }

  double toReturnS = (180/PI)*sqrt(toReturn);
  return toReturnS;
}

// shortes distance between two angles
double Molecule::getDistanceBetweenAngles(double &a1, double &a2)
{
  double angle = fmod(abs(a1-a2), 2.0*PI);

  if (angle > PI) angle = 2.0*PI - angle;

  return angle;
}